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Abstract 

We study strongly attractive fermions in an optical lattice superimposed by a trapping potential. We 
calculate the densities of fermions and condensed bound molecules at zero temperature. There is a 
competition between dissociated fermions and molecules leading to a reduction of the density of fermions 
at the trap center. 
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1. Introduction 

Ultra-cold gases present many-body systems with a remarkable tunability of its parameters, providing 
us with a platform for the investigation of numerous properties of complex quantum systems. For instance, 
it has become possible to tune the interaction strength between atoms over a wide range by Feshbach 
resonance The investigation of ultra-cold Fermi gases started shortly after the discovery of Bose- 
Einstcin condensation of bosonic atoms Q when quantum degeneracy in a gas of fermionic atoms was 
obtained Q . For fermions with weak attraction we have the celebrated phenomenon of Cooper pairing Q . 
By tuning to strong attraction one enters the regime of diatomic molecules, whose size is much smaller 
than that of Cooper pairs. These molecules then may condense into a Bose-Einstein condensate of a 
hard-core Bose gas. The crossover from the weakly interacting BCS regime to the strongly interacting 
BEC regime has been the subject of theoretical studies in 0, Q. The bosonic molecules, formed of pairs 
of fermionic atoms, were produced experimentally in a trapped system as well as in an optical lattice 
i- 

Most investigations were based on continuous Fermi gases. The more recently introduced optical 
lattices in ultra-cold gases may have a number of interesting effects on the Fermi gas. First of all, 
the dispersion of the atoms will be changed by the lattice. Moreover, the interaction between the atoms 
has a strong effect on the quantum states of the Fermi gas by allowing, for instance, to form Mott states 
[loj . This could mean that the molecular gas in the BEC regime does not condense but becomes a Mott 
state. In the BCS regime, on the other hand, the effect of an optical lattice is not so dramatic because 
the Cooper pair radius is much larger than the lattice spacing. Consequently, the BEC-BCS crossover 
can be much richer in the presence of an optical lattice. 

In this paper, we study strongly attractive fermions in an optical lattice superimposed by a trapping 
potential. We first show that the phase diagram of tightly bound fermions contains a Bose-Einstein 
condensed phase and a Mott insulating phase of such molecules. Then we study the system in a harmonic 
trap and calculate the density of unpaired fermions in the presence of a condensate state of molecules. 
We show that there is a competition between paired fermions and unpaired fermions which leads to a 
reduction of the density of unpaired fermions at the center of the trap. 
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2. Model 



The simplest lattice model of an interacting Fermi gas is the Hubbard model It describes the 
competition between the kinetic energy of the fcrmions and a local interaction and provides a phase 
diagram that includes a Fermi liquid and a Mott insulator. Originally introduced for a repulsive Fermi 
gas (e.g. electrons in a metal), the model can also be used for neutral fermionic atoms with attractive 
local interaction. For very strong attraction, however, the local interaction is insufficient to describe the 
physics of the Fermi gas. This due to the fact that strong attraction causes pairing of fcrmions to local 
molecules. On the other hand, the kinetic term in the Hubbard model allows only individual tunneling 
of fermions. This means that the tightly bound molecules must dissociate into independent fcrmions 
in order to tunnel in the optical lattice. The associated energy of such a process is of the order of the 
attractive interaction. (Actually, the effective tunneling rate is ~ It} jU @, where t is the tunneling rate 
of individual fermions and \J is the strength of the local attraction.) This means that the motion of 
molecules is strongly suppressed in the Hubbard model [l^ . On the other hand, there is no reason for 
the molecules not tunnel freely in the optical lattice because they only have to obey the Pauli principle. 
The solution of this problem is an extension of the attractive fermionic Hubbard model that includes an 



additional kinetic term for the bosonic molecules |13| - |16| |. The Hamiltonian for an attractive Fermi gas 



in a d-dimensional optical lattice is a molecular fermionic Hubbard (MFH) model and reads 

o-=t4 {r,r') {r.r') r o-=t4 

Here, Cr^a (cJo-) is the annihilation (creation) operator for particles at lattice site r. The index a =t,i 
represents two hyperfine states of fermionic atoms, e.g., or ^Li. In this work we consider only the 
symmetric case where the number of fermions in each component is the same. Nearest-neighbor tunneling 
of the individual fermions is described by the parameter t. There is also a term with parameter J which 
is understood as a tunneling term of dressed fermionic pairs [H, 15|. U ^ Ubg — -idg^/S accounts for an 



effective local attractive interaction between fermions with the detuning S, coupling between fermions 
and molecules g and Ubg = ^Trll^ab/m with the background scattering length Ub and m is the mass of 
fermions. In the strong coupling regime \S\ ^ g and ab is small and positive [l^. So that for strong 
attractions U < J. A parabolic trapping potential Vr = 7(2;^ + j/^) (we will study a 2D system) can be 
combined with the chemical potential fi^ to fXra- = fJ-a ~ Vr, which controls the number of particles in a 
grand-canonical ensemble, 7 is the strength of the trap. 



3. Functional integral representation 

The grand-canonical ensemble is given by the partition function as a functional integral with respect 
to the Grassmann fields, Z = J e'^'^'^^'^^Dli), with the action [l^ 



dT 



(2) 



where e should be sent to +0 at the end of calculations, r is the imaginary time, the symbol dripr denotes 
the formal lime_j._|_o('0T-i-£ — V'r), P = {kBT)~^ is the inverse temperature. H f{ipT+e,'4'T) is obtained 
by replacing creation and annihilation operators in the Hamiltonian Hf hy Grassmann fields -0 and 
respectively. We will perform a Hubbard- Stratonovich transformation to decouple the fourth order terms 
at the expense of introducing new complex fields. We use the identity 
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Figure 1: Phase diagram for t = 0, U = and fcsT = 0. The three phases are the BEC with a non-vanishing condensed 
density, the MI states with vanishing condensed density and with one molecules per lattice site and the empty phase, fi 
and J are in arbitrary energy units. 



(3) 



where/ = 6{J—\U\),Vr,r' ~ jj_S\r-r'\.a — J<>r,r' + i^+I)USr,r'- The subsequent integration over Grassmann 
fields leads to the effective action in terms of the two complex fields (f), x ^-nd their conjugate fields 0, x 



dT 
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with an inverse Nambu-Gor'kov propagator 



dr — l-Lf - t 20 + Ix 



(4) 



(5) 



Here i(f) + /x is a diagonal matrix with i(j)rr + IXrr on its diagonal, fia- is also a diagonal matrix with 
elements jira- i is a hopping matrix with elements t/2d. 

The total density of fermions can be calculated as a derivative of the free energy F = \nZ/(3 with 
respect to the chemical potential, n^o- = dF/dfira- In the strongly interacting regime in an optical lattice 
we expect that the size of a molecule is within one well of the optical lattice and the internal structure is 
irrelevant. The corresponding creation operator is cj^^c^^j^. Then the molecules can be treated as point-like 
objects. The scattering length of fermions is assumed to be less than the lattice spacing of the optical 
lattice, since the size of the molecules is of the order of the scattering length. 

In order to calculate the condensate density we consider the presence of "off-diagonal long-range 
order" [H 

no = , lim - / dr hm (?/irr-|-etV'r'rtV'rT+£i'0r'ri)- (6) 



\r—r' I - 



/3 



>+0' 



limi, 



It can be shown that the latter is related to the correlations of the complex fields, tiq 

where (...) =Z-i/---e"^="-D[0,x] El- 

In the limiting case of strong attraction all fermions are paired and the tunneling of molecules occurs 
without their dissociation. In this case U can be neglected and the atomic tunneling rate t is not effective. 
Without a harmonic trap (7 = 0), the integration over the complex molecular fields 0, x can be performed 
in saddle-point approximation [20|. The phase diagram is depicted in Fig. [TJ There are three phases: 
the BEC of molecules with the condensed fraction no = ( — 4/i^)/(4J^), the MI state with one particle 
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per site and no = and the empty phase (cf. 116|). Fluctuations around the saddle point provide the 
low-energy excitations of the bosonic molecules. They are gapless in the BEC phase 

Efc = ^Ji.Pgkna+4glfi^, (7) 

where gk — 1 — 1/ d Y^i^i cos(fcia) is the dispersion of the free Bose gas. On the other hand, the excitations 
of the MI state have a gap A = 2^ — J > 

Ek^A + Jgk. (8) 



4. Ground state of the trapped system 

The approximate ground state of the system is obtained as the saddle point of the action in Eq. Q 
with respect to the fields (i.e., by solving the equation 6Seff ~ 0). Here we assume that the solution is 
static (i.e. independent of r). Moreover, we also use the fact that a slowly varying field (p in space of a 
trapped condensate can be approximated as 

V~^,4>r' ~ i>(j)r + Jb^ ^ i5rr' - Xr' ) <l>r' , (9) 



where b ^ = (1 + I)U, J is a tunneling matrix with elements l/2d. Then we obtain from SSeff = the 
following equations for the fields (j), 4> 

Jb'^ - - - - 1 

- ^(0r+ae, " 2(^r + (f>r-aei) = -bcj)r " '-^rr,ll{'-^n) , (10) 



1=1 

- ^('/>r+ae, ^ 20^ + 0r-ae, ) = + i^ ^ G„.,22 (l^ri ) , (H) 

1 = 1 ^ n 

where ei is the Cartesian lattice unit vector in direction I. Grr' ,11-12-21-22 are defined as the 2x2 block 
structure of 

G=f ^-'12 V (12) 

\ ^rr',21 ^rr'.,22 J 

r ~ 1,2,..., N'^, where G is the inverse matrix of 

G-'^( zo; + t y 

\ lUn - l-Lf l<p + Ix J 

= 7r(2n + 1)113 is a Matsubara frequency originating from the Fourier transformation dr — > — iw„. 
The equations of the other complex field x read 

Xr = -/t/-^^G„,ii(w„), (14) 

n 

Xr = lU^^^ Grr,22{^n) ■ (15) 
n 

Eqs. (fTO]) and (fTTj) are analogous to the Gross-Pitaevskii equation for a Bose gas and provide a macroscopic 
wave function of the condensate molecules. The left-hand side of Eqs. (fTO)) and pT|) can be understood 
as a lattice Laplacian acting on the field 

d 

S/'^cj)^ = ^(0r+aei - 2(f>r + 4>r-aei) ■ (16) 
1=1 
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It should also be noticed that the additional condition V^(/) = gives the Thomas-Fermi approximation. 
Using the effective action in Eq. (|4]) and differentiating it with respect to the chemical potential we get 
the total densities ^ ^ 

n-r^ = ^^^Grr.l2{(^n), rirl = - Grr.2l{i^n) ■ (17) 

n n 

In mean- field approximation (cf. Ref. (20l |) this gives for the densities of unpaired fcrmions 

4^ = n,, - n.,n., - no,., 4^ = n., - n.,n., - no,. (18) 

Here Uf measures the presence of a fermion and absence of a pair of fcrmions. rir^nri gives the probability 
to find two fcrmions at a lattice site r. Finally, tiq, is the product of anomalous averages ~ {ip-i-ipi) (V'iV't) 
and thus can be associated with the condensed density 

no, = 4 X! Grr,ll{uJn)^ ^ G„.,22(w,i) \i(j>r + Xr? ■ (19) 
n n 



5. Numerical calculation 



We calculate numerically the densities of Fermi gas in a 2D dimensional N x N (N ~ 30) sites in an 
optical lattice. We explored the case J > \U\ that corresponds to the strongly interacting regime [isj . 
We choose the lattice constant a = 0.3 in our calculation. Equations (jlOp and (fTT|) are elliptical equations 
and similar to the Poisson equation. A well known technique to solve the Poisson equation numerically 
is the relaxation method, where an artificial time dependent term is added [21|. In two dimensions, for 
equations (|10I11I14I15P we obtain the following relaxation equations 



dt 



dt 



2d - 2d Xi j 

+ + - + 4>i.j-i - —(t>i,j - ii^iTf, (20) 



2d 



m lU 

. 2d Xi,3 
'j62 lU 



dXr 

dt 

dXr 

dt 



= -Xr 



lU 



rr,22 



(l^n) 



(21) 
(22) 

(23) 



,2,...,N is a Cartesian mapping to r = 1,2,...,N'^ that describes the coordinates of the 
Eqs. (|20| and (|2T|) arc evolved through a splitting operator technique 



2lj . The corre- 



where i,j = 1, 
square lattice. 

spending terms are evolved by a Crank-Nicolson algorithm while the remaining ordinary equations 
are evolved by Euler method. The equations are evolved until 0, x,x no longer changes with time, 
i.e., d(j)/dt = d(j>/dt = dx/dt = dx/dt = and thus the set of equations are equation is satisfied. Equa- 
tions ()20|21|22|23p where discretized with time steps ranging from 0.01 to 0.001. We start the evolution 
of eqs. ([20|) and ((2T|) from an initial Gaussian ansatz for (f>, (p and X = X = i- We then apply a nu- 
merical diagonalization to equations described in the appendix [Appendix A| and obtain the summations 
1/(3 J2n G'rr,ii(w„) and 1/(3 J^n Grr.ui^n)- Thcu wc cvolvc <f>, (j) with the V'^ terms of eqs. (|20l2ip . evolve 
X, X in cq. (|22l23p and finally evolve <f>, (j> with the remaining terms of cq. (j20l2ip . These updated values 
will then be used as new ansatz for numerical diagonalization and the circle is repeated. Convergence 
was achieved after t ^ 10, or 1000 steps. The above numerical procedure gives the complex fields (f> and 
X as a function at each site r, i.e., the field at each point in the Cartesian grid i, j in the plane x, y, where 
x = {i-N/2- l/2)a, y=U- N/2 - l/2)a, i,j = 1, 2, ...N. 

First we consider a situation when the Hamiltonian contains only the tunneling of paired fermions, 
i.e., for t = Q and U = Q. This corresponds to the limiting case of the strongly interacting regime [3, llBl- 
We also fix 7 = 0.05J and large chemical potential fi^ = J (large number of fermions). In Fig. [2] a Mott 
plateau of the paired molecules is shown. The condensed molecules (no) form a superfluid shell, which 
assembles around the Mott plateau (cf. [22j). 
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Figure 2: Mott state in the trap for values of f/J = and (7 = 0. The chemical potential is fi^ = J- The superfluid shell 
assembles around the Mott plateau. The density profiles are along the lines y = a/2 and x = a/2. 

Now we fix the smaller chemical potential fig- ^ and nonvanishing t and U. In Fig. [3] the condensed 
density as well as the density of dissociated atoms are plotted for U = 0.5 J. We notice that there is a dip 
at the center of the trap for the dissociated atoms. If the condition for the strongly interacting regime 
is not satisfied (e.g., U = 1.5 J in Fig. 2]), then this dip is less pronounced. In Fig. [5] we also plot the 
densities in the Thomas-Fermi approximation (by putting terms on the left hand sides in Eqs. pOI) and 
((TT|) to zero) that does not reveal the dip formations as compared to Fig. [3] 

Thus, the contribution of the nonlocal nabla terms in Eqs. (|10p and (jlip . which originate from the 
tunneling of the molecules in Eq. ([T]), is essential for the appearance of the dip. The tunneling of such 
dressed molecules is due to an optical lattice and their fermionic nature [Tsj . So the appearance of the dip 
can be explained tentatively by the Pauli principle which acts between paired and unpaired fermions and 
causes the paired fermions to repel the dissociated ones (similar mechanism acts for unbalanced Fermi 
mixtures |23[). Eq. ([T]) describes a one band fermionic model and thus is valid for a deep optical periodic 
potential. For a shallow lattice the Pauli principle is less efficient and we expect the dip formation to 
be less pronounced. The Pauli principle is also less efficient if the density of fermions is small. We may 
expect that in the dilute regime the role of the lattice is almost irrelevant. Indeed, if we decrease the 
chemical potential, i.e., by going into the dilute regime, the dip softens (cf. Fig. [5]), while increasing the 
chemical potential makes the dip more profound (cf. Fig. [7]). Thus, the role of the lattice is crucial in 
the formation of the dip in the density of unpaired fermions at the trap center. 

In conclusion, we have studied strongly attractive fermions in a 2D optical lattice superimposed by 
a trapping potential. We solved numerically the field equations obtaining non-homogeneous densities of 
condensed molecules and dissociated fermions. The latter reveal an effect of dip formation at the trap 
center. This can possibly be observed in future experiments and may serve as a signature of approaching 
of the strongly coupling BEC regime of the BEC-BCS crossover in a lattice. 

O.F. and K.Z. acknowledge support from DAAD. A.G. thanks support from CAPES (Brazil) and 
FAPESP/CNPq (Brazil). 

Appendix A. Summation over Matsubara frequencies 

Here we show how we perform the summation over Matsubara frequencies. We transform iuJn z. 
G~^{z)G~^{z) is a Hermitian matrix (iA = 10 + x = + is real, x is complex): 

^ ^''^ - [ UA - lAi ^ A2 + z2 _ ^2 _ £^ _ ^£ _ I • (A-l) 
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Figure 3: The condensate density and the density of dissociated atoms for t/ J = 0.3, U = 0.5 J. There is a dip at the center 
of the trap in the case of the dissociated atoms. It maybe explained by the Pauli exclusion principle acting between paired 
and unpaired fermions. The density profiles are along the lines y = a/2 and x = a/2. 



o 




Figure 4: The condensate density and the density of dissociated atoms for t/J = 0.5, U = 1.5 J. There is no dip at the 
center of the trap in the case of the dissociated atoms. The density profiles are along the lines y = a/2 and x = a/2. 



Any Hermitian matrix can be diagonalized by a unitary matrix: 

G-^{z)G-\z)^U'^\{z)U, U^ = tj-\ 
where A is a diagonal (real) eigenvalue matrix. The equation for A of Eq. (IA.1[) reads 



Reversing Eq. (|A.2p we get 



\ 2 2 

Afc = z - zi.. 



G{z) = il^— ^UG~^{z). 

Z - Zi, 



(A.2) 

(A.3) 
(A.4) 
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Figure 5: The condensate density and the density of dissociated atoms in Thomas-Fermi approximation for the same values 
of tj J and U as in Fig. [S] The Thomas- Fermi approximation can not reproduce the dip formation. The density profiles 
are along the lines y = a/2 and x = a/2. 
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Figure 6: The condensate density and the density of dissociated atoms for the same values of t/J and U as in Fig. [3] but 
for lower densities with fia = — O.IJ. The dip disappears since at lower densities the role of the lattice is less profound as 
for higher densities. The density profiles arc along the lines y = a/2 and x = a/2. 



For the well behaved g{x) (see [17| ) 

m2 -X- r? 



1 



1 
2TTi 



dz 



9{z)f{z) ^ g{zk)f{zk) - g{-Zk)f{-Zk) 
- zl 2zk 



(A.5) 



where f{z) = l/{e^^ + 1). This works for the blocks 11 and 22 of the Green matrix in Eq. (O, since the 
function under the integral is well-behaved then. For the blocks 12 and 21 the sum is formally divergent 
since for large w„ it behaves as ^ 1/wn. To cure the problem we must introduce convergent factors 
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Figure 7: The condensate density and tlie density of dissociated atoms for the same values of f/J and U as in Fig. [3] but 
for higher densities with ficr = O.IJ. The dip deepens since at higher densities the role of the lattice is more profound. The 
density profiles are along the lines y = a/2 and x = a/2. 



Using the above formula and the remark we get 

i ^ GrrM = ^^^f^{tanh[-/3zfe/2]G,-U0) - Zk^a2}■ (A.6) 



References 

[1] E. Timmcrmans, P. Tommasini, M. Hussein, A. Kerman, Phys. Rep. 315, 199 (1999); I. Bloch, J. 
Dalibard, W. Zwerger, Rev. Mod. Phys. 80, 885 (2008). 

[2] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, E. A. Cornell, Science 269, 198 
(1995); K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, W. 
Kctterlc, Phys. Rev. Lett. 75, 3969 (1995). 

[3] B. DeMarco, D. S. Jin, Science 185, 1703 (1999). 

[4] J. Bardeen, L.N. Cooper, J.R. Schrieffer, Phys. Rev. 108, 1175 (1957). 

[5] A.J. Leggett, J.Pliys. (Paris), CoUoq. 41, 7 (1980); C.A.R. Sa de Melo, M. Randeria, J. R. Engelbrecht, 
Phys. Rev. Lett. 71, 3202 (1993); J. R. Engelbrecht, M. Randeria, C.A.R. Sa de Melo, Phys. Rev. B 
55, 15153 (1997); M. Holland, S. J. J. M. F. Kokkelmans, M. L. Chiofalo, R. Walser, Phys. Rev. Lett. 
87, 120406 (2001); Y. Ohashi, A. Griffin, Phys. Rev. Lett. 89, 130402 (2002); Y. Ohashi, A. Griffin, 
Phys. Rev. A 67, 033603 (2003). 

[6] P.Nozicres, S.Schmitt-Rink, J. Low Temp. Phys. 59, 195 (1985). 

[7] C. A. Regal, Ch. Ticknor, J. L. Bohn, D. S. Jin, Nature 424, 47 (2003); I. D. Prodan, M. Pichler, M. 
Junker, R. G. Hulet, Phys. Rev. Lett. 91, 080402 (2003); M. Grcincr, C. A. Regal, D. S. Jin, Nature 
426, 537 (2003); S. Jochim, M. Bartenstein, A. Ahmeyer, G. Hendl, S. Riedl, C. Chin, J. Hecker 
Denschlag, R. Grimm, Science 302, 2101 (2003). 

[8] T. Stoferle, H. Moritz, K. Giinter, M. Kohl, T. Esslinger, Phys. Rev. Lett. 96, 030401 (2006); J. K. 
Chin, D. E. Miller, Y. Liu, C. Stan, W. Setiawan, C. Sanner, K. Xu, W. Ketterle, Nature 443, 961 
(2006). 



9 



[9] M. Greiner, I. Bloch, O. Mandel, T. W. Hansch, and T. Esslinger, Phys. Rev. Lett. 87, 160405 (2001); 
W. Ketterle, M. Inguscio, C. Salomon, Ultra-cold Fermi Gases, Proceedings of the Varenna "Enrico 
Fermi" Summer School (2007); S. Giorgini, L. P. Pitaevskii, S. Stringari, Rev. Mod. Phys. 80, 1215 
(2008). 

[10] M. P. A. Fisher, P. B. Weichman, G. Grinstein, D. S. Fisher, Phys. Rev. B 40, 546 (1989); D. Jaksch, 
C. Bruder, J. I. Cirac, C. W. Gardiner, P. ZoUer, Phys. Rev. Lett 81, 3108 (1998); M. Greiner, O. 
Mandel, T. Esslinger, T. W. Hansch, I. Bloch, Nature 415, 39 (2002). 

[11] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, U. Sen, Adv. Physics bf 56, 243 (2007). 

[12] G. Orso, L.P. Pitaevskii, S. Stringari, M. Wouters, Phys. Rev. Lett. 95, 060402 (2005). 

[13] K. Ziegler, Laser Physics 15, 650 (2005). 

[14] Y. Ohashi, Phys. Rev. A 78, 063617 (2008). 

[15] L.-M. Duan, Phys. Rev. Lett. 95, 243202 (2005). 

[16] L. D. Carr, M. J. Holland, Phys. Rev. A 72, 031604(R) (2005). 

[17] A. Altland, B. Simons, Condensed Matter Field Theory (Cambridge University Press, New York, 
2006). 

[18] C. N. Yang, Rev. Mod. Phys. 34, 694 (1962); A. J. Leggett, Quantum Liquids (Cambridge University 
Press, New York, 2006). 

[19] Ch. Moseley, O. Fialko, K. Ziegler, Ann. Phys. (Berlin) 17, 48 (2008). 

[20] O. Fialko, Ch. Moseley, K. Ziegler, Phys. Rev. A 75, 053616 (2007); O. Fialko, K. Ziegler, J. Phys. 
B 43, 065304 (2010). 

[21] W. H. Press et al.. Numerical Recipes, 2nd ed. (1992) 

[22] Y. Fujihara, A. Koga, N. Kawakami, Physica B 404, 3324 (2009). 

[23] M. Iskin, C. A. R. Sa de Melo, Phys. Rev. A 77, 013625 (2008); M. Iskin, C. J. Wilhams, Phys. Rev. 
A 78, 011603(R) (2008). 

[24] R.B. Diener, R. Sensarma, and M. Randeria, Phys. Rev. A 77, 023626 (2008). 



10 



